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Abstract 


Jointly  Constrained  Bilinear  Programming: 
The  Linear  Complementarity  Problem 


Faiz  A.  Al-Khayyal 

School  of  Industrial  and  Systems  Engineering 
Georgia  Institute  of  Technology 
Atlanta,  Georgia  30332-0205 


-We- investigates; ref inements  to  an  existing  nonconvex  programming 
algorithm  that  exploit  the  special  structure  of  linear  complementarity 
problems.  Wf  prove/. that  the  working  bases  in  the  linear  programming 
subproblems  can  be  reduced  from  3nx3n  to  nxn.  In  addition,  we-  show/^that 
the  procedure  (in  general,  infinitely  convergent)  is  finite  under  a 
nondegeneracy  assumption.  The  procedure  compares  favorably  with  two 
recently  proposed  algorithms  and  is  competitive  with  a  third  related 
method.  ^ 
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Jointly  Constrained  Bilinear  Programming: 

The  Linear  Complementarity  Problem 

Faiz  A.  Al-Khayyal 

School  of  Industrial  and  Systems  Engineering 
Georgia  Institute  of  Technology 
Atlanta,  Georgia  30332-0205 

In  a  recent  paper,  Al-Khayyal  and  Falk  [1]  describe  a  branch-and- 
bound  algorithm  for  finding  a  global  solution  to  the  nonconvex  program 

T 

minimize  f(x)  +  x  y  +  g(y) 

subject  to  (x,y)  e  S  n  £2  (1) 

where  f  and  g  are  convex  over  the  feasible  region,  which  is  defined  by 
the  intersection  of  a  nonempty,  closed,  convex  set  S  and  a  compact  hyper- 
rectangle  0.  For  an  appropriate  choice  of  0  the  above  problem  contains 
as  a  special  case  the  well-known  linear  complementarity  problem  (LCP)  of 
finding  a  real  n-vector  x  such  that 

Mx  +  q  >  0,  x  >  0,  x^(Mx+q)  =  0,  (2) 

where  M  is  a  given  real  square  matrix  and  q  is  a  given  real  n-vector. 

The  relationship  of  (2)  to  (1)  is  apparent  from  the  following  indefinite 
quadratic  programming  formulation  of  (2): 


minimize 


T 

x  y 


subject  to  -Mx+y  =  q,  x  >  0,  y  >  0 


(3) 


To  solve  (3)  by  the  branch-and-bound  method  in  [1]  we  need  a  compact 
hyperrectangle  ft  that  contains  at  least  one  complementary  solution.  Such 
a  set  is  easy  to  construct  in  this  case  because,  when  the  LCP  has  a  solu¬ 
tion,  at  least  one  complementary  solution  is  a  basic  feasible  solution  of 
the  system  -Mx+y  =  q.  Hence  an  appropriate  ft  may  be  constructed  by 
bounding  all  basic  solutions  using  the  following  result  [6,  Lemma  2.1]. 


Lemma  1 


Let  x  be  a  basic  solution  of  the  mxn  system  Ax  =  b.  Then 


|xj  |  <  m!  o”-1 0 


where 


max { | a  J }  and  g  »  max { | b . | } . 

i,j  J  j  2 


Tighter  bounds  can  be  obtained,  in  the  case  when  the  set  D  =  {x:  Mx+q  > 


0,  x  >  0}  is  bounded,  by  solving  the  n  linear  programs  max{Xj :  xeD}  for 


j  =  l,...,n.  Once  upper  bounds  on  x  are  determined,  similar  bounds  for 
y  =  Mx+q  are  easily  computed. 

In  this  paper  we  investigate  refinements  to  the  method  in  [1]  that 
exploit  the  structure  of  problem  (3).  In  particular,  we  show  that  the 
size  of  the  working  bases  in  the  linear  programming  subproblems  can  be 
reduced  from  3nx3n  to  nxri,  and  that  convergence  is  finite  for  linear 
complementarity  problems  with  nondegenerate  complementary  solutions.  We 
also  compare  the  procedure  to  three  recent  algorithms.  The  preliminary 
findings  are  that  the  proposed  procedure  is  faster  than  the  cutting  plane 


*» 


3 


method  of  Ramarao  and  Shetty  [7],  and  can  solve  more  problems  than  the 
piecewise  linear  equation  approach  of  Solow  and  Sengupta  [8].  Only  the 
enumeration  procedure  of  Al-Khayyal  [3]  is  potentially  faster,  but  the 
two  methods  appear  to  give  comparable  execution  times  on  some  problems. 
For  brevity,  we  assume  knowledge  of  the  algorithm  in  [1]  and  will  make 
reference  only  to  components  of  the  latter  procedure  that  are  provably 
refined  here. 


1.  Algorithmic  Refinements 

For  the  general  problem,  the  algorithm  in  [1]  branches  into  four 

subproblems  at  each  stage.  These  problems  are  defined  by  partitioning 

the  hyperrectangle  SI  into  four  subsets  in  the  following  way.  Given  a 

point  (x,y)efl  =  {(x,y):  l  <x<u  ,  l  <y<u  },  we  choose  an  index  I  to  speci- 

x  x  y  y 

1  4 

fy  the  partition  of  SI  into  the  four  sets  51 . SI  ,  where 


n1  «  ((x,y):  *Xi<xx<v  ^y^ylt 


V<Xi<Ux  *  V<yi<Uy  » 


*i  yi 


2  —  — 

SI  =  {(x,y):  xI<xI<ux  ,  l 


I  yI 


l  «  <u  ,  l  *y. <u  ,  i*I } , 
x,  i  x  ’  y_.  •'i  y  *  J’ 


*i  yi 


A  =  {(x,y):  xI<xI<ux  ,  y^y^u^  , 

l  <x  <u  ,  l  <y. <u  ,  i*l}, 
Xi  1  Xi  yi  1  yi 


,  v.V  ’.Wv’' 


£2  =*  {(x,y):  l x  OCjttj,  yI<yI<uy^» 

i  <x  <u  ,  l  <y. <u  ,  i*l}. 
xi  1  xi  yi  1  yl 


Each  of  Che  four  subsets  is  a  candidate  for  future  partition  in  the 
search  for  an  optimal  solution.  Without  loss  of  generality,  we  may 
assume  that  u  >  l  >0  and  u  >  l  >  0. 

xx  y  y 

A  solution  to  the  linear  complementarity  problem  (2)  must  satisfy 

Xiyi  *  ®  f°r  **  Hence»  the  index  I  is  chosen  such  that  x^y^.  >  0  at  a 

feasible  noncomplementary  solution  (x,y).  The  set  £2  chosen  for  partition 

is  the  one  among  all  candidate  subsets  that  produces  Che  lowest  value  of 

xTy,  where  each  subset  has  a  point  (x,y)  associated  with  it.  Clearly, 

3 

the  subset  £2  can  be  eliminated  from  further  search  (fathomed)  since  it 

cannot  contain  a  complementary  solution.  We  show  in  the  next  section 

that  the  set  £2  can  be  fathomed  if  either  l  >  0  or  l  >0  together  with 

XI  yI 

x  y  >  0.  Thus,  I  is  a  partitioning  index  for  £2  only  if  2  -  i  -0. 

XI  yI 

In  this  way,  the  active  nodes  (corresponding  to  subsets  that  potentially 
contain  complementary  solutions)  are  guaranteed  to  satisfy  l  l  =0  for 


Xi  yi 


all  i. 


The  greatest  computational  savings  are  realized  from  the  refinement 
of  the  subproblem  solved  at  each  node.  Let  £2  be  the  active  partition 
associated  with  an  arbitrary  node.  The  subproblera  solved  at  that  node  is 


defined  as 


linear  program  involving  3n  decision  variables  (2n  are  lower  and  upper 
bounded,  and  n  are  unrestricted)  and  3n  constraints  (n  are  equality  and 


and  N  =*  {l,2,...,n}.  When  S  nc  *  0,  let  (x,y)  denote  an  optimal  solution 
to  problem  (5)  and  let  v  =  ^(x,y).  Notice  that  C=ft  and  that  i))(x,y)  > 

4»(x,y)  for  all  (x,y).  In  addition,  \Ji(x,y)  >  0  for  (x,y)  >  (0,0).  We 

show  the  equivalence  of  problems  (4)  and  (5)  for  active  partitions  ft,  but 
first  we  prove  the  following  simple  Lemma. 

Lemma  2  If  S n  ft  contains  a  complementary  solution,  then  v  =  v  =  0. 

*  * 

Proof  Suppose  (x  ,y  )eSnft  is  a  complementary  solution.  By  con- 

T 

struction,  we  have  l  >  0  and  l  >0  such  that  l  i  *  0.  It 

x  y  x  y 

k  k 

follows  that  (x  ,y  )eC.  Because  ip  is  the  convex  envelope  of 
T 

x  y  over  ft,  we  have 

0  *  (x  )Ty  >  <J)(x  ,y  )  >  <(>(x,y)  >  i|>(x,y)  >  i^^.y)  >0.  B 

Problems  (4)  and  (5)  are  equivalent  in  the  sense  that  either 
the  optimal  solution  of  one  solves  the  other,  or  S n ft  does 
not  contain  a  complementary  solution. 

We  first  show  that  either  the  optimal  solutions  of  problem 
(4)  solve  problem  (5)  or  S  n ft  does  not  contain  a  complemen¬ 
tary  solution.  If  (x,y)eC,  then  $(x,y)  =  \Ji(x,y).  Hence, 

v  *  min{$(x,y):  (x,y)eS  nc) 

>  min{iji(x,y) :  (x,y)eSn  C) 

>  min{4i(x,y):  (x,y)eSnft|  =»  v 

and  therefore  (x,y)  solves  problem  (5).  Now  suppose  that 


Theorem  3 


Proof 


(x,y)^C.  Then  v  >  0  and,  by  Lemma  2,  S  n ft  does  not  contain 
complementary  solution. 

To  complete  the  proof,  we  now  show  that  either  the  opti 
mal  solutions  of  problem  (5)  solve  problem  (4)  or  S  n ft  does 
not  contain  a  complementary  solution.  The  latter  holds  when 
S  nc  =  0,  since  every  complementary  solution  in  S  n  jj  must 
necessarily  be  in  C.  Now  assume  that  S  nc  *  0.  We  consider 
two  cases:  when  v  =  0  and  when  v  >  0.  First  assume  that  v 
0  and  recall  that  (x,y)  is  feasible  to  problem  (4).  We  have 

0  =  v  =  $(x,y)  >  (j>(x,y)  =  v  >  0. 

Hence  <J> ( x , y)  =  0  and  (x,y)  solves  problem  (4).  Finally, 
assume  that  v  >  0.  We  now  prove  that  v  >  0  and  invoke  Lemma 
2  to  claim  that  S  n  ft  does  not  contain  a  complementary  solu- 
tion.  By  contradiction,  assume  that  v  >  0  and  v  =  0.  Then 

0  »  $(x,y)  >  ^(x,y)  >  0 

and  (x,y)  is  feasible  to  problem  (5);  otherwise,  (x,y)£C  and 
<j>(x,y)>0.  Hence, 

0  =  ij/(x,y)  >  ijKx.y)  ■  v  >  0 

which  is  a  contradiction.  Thus  we  conclude  that  v  >  0.  I 


By  exploiting  the  foregoing  theorem,  the  branch-and-bound  algorithm 
[1]  involves  linear  programming  subproblems  with  working  bases  of  order 
n.  This  is  achieved  by  scaling  the  variables  in  problem  (5).  Specifi¬ 
cally,  for  each  ieN,  let 


\\hf 


and  y! 


/u 

-l  \ 

X,  1 

I  _i 

i  1 

1  u 

u  I 

\  Xi 

In  terms  of  the  new  variables,  problem  (5)  becomes 

T  T 

min  c  x'  +  d  y* 


subject  to 


Ax'  +  By'  =  q 
x’  +  y'  <  b 


y'  >  *y. 


where  the  parameters  c,d,A,B,b,£  ,  and  l  ,  are  easily  derived  from  the 

x  y 

coordinate  transformation  (6).  Problem  (7)  need  not  be  solved  to  comple¬ 
tion  if  a  complementary  solution  is  uncovered.  It  is  easy  to  recognize 

T 

when  such  a  solution  is  encountered  because  x  y  =  0  if  and  only  if 
( x ' ) Ty '  =  0. 

The  linear  program  (7)  has  n  equality  costraints,  n  generalized 
upper  bounding  (GUB)  constraints,  and  2n  lower  bound  (LB)  restrictions  on 
the  decision  variables.  Recall  chat  both  GUB  and  LB  constraints  can  be 


handled  separately  (see,  e.g.,  Chapter  6  in  Lasdon  [4]).  Therefore,  all 


pivot  computations  can  be  performed  using  only  an  nxn  basis  matrix,  which 
is  a  considerable  improvement  on  the  3nx3n  basis  that  arises  in  the  gene¬ 
ral  case. 


2.  Finite  Convergence 

For  the  general  problem  (1),  the  branch-and-bound  algorithm  [1]  can 
only  be  guaranteed  to  converge  in  the  limit.  In  this  section  we  prove 
that  the  process  is  finite  for  linear  complementarity  problems  whose 
solutions  are  nondegenerate  extreme  points  of  S  njj.  We  first  show  that 
the  procedure  always  branches  into  three  subproblems  and  then  argue  that 
all  paths  in  the  branch-and-bound  tree  are  fathomed  after  a  finite  number 
of  branches  under  the  nondegeneracy  assumption. 

Let  ft  represent  the  initial  hyperrect angle  and  let  C  be  the  associ¬ 
ated  set  in  problem  (5).  Thus,  £  *  £  =  0  and  i|;(x,y)  =*  0  for  all  (x,y). 

x  y 

Suppose  (x,y)eS  nC  such  that  x^y^.  >  0*  Note  that  (x,y)  solves  problem 
1  4 

(5).  Let  ft  ,...,ft  be  the  partition  of  ft  defined  above,  and  let 
14  14 

C  ,...,C  and  ip  , . . . ,  ip  be  the  associated  sets  and  objective  functions, 

respectively,  in  problem  (5).  For  j  =  1,...,4,  let  (r'.r1)  and  (uJ,uJ) 

x  y  x  y 

denote  the  lower  bound  and  upper  bound  vectors  in  ft^ . 

3  3  3 

Since  (£  ,£  )  >  (0,0),  the  set  ft  cannot  contain  a  complementary 

XI  yi 

solution;  therefore,  we  fathom  that  subset.  For  j  =  1,2,4,  if  S  n  =  0 
then  by  Theorem  3  the  set  ft^  also  does  not  contain  a  complementary  solu¬ 
tion  and  is  consequently  fathomed.  Otherwise,  let  (x'^.y'')  denote  an 
optimal  solution  to  problem  (5)  with  ip  =  ^  and  C  =  G^.  For  conveni¬ 
ence,  we  shall  henceforth  designate  such  problems  as  (5-j).  Now,  for  j  = 
2,4,  if  ijp(xj  .yJ)  >  0  then  the  subset  ft^  is  fathomed  by  Lemma  2.  Hence, 


set  ft^  is  a  candidate  for  further  partition  only  if  i|^(x^,y^)  *  0.  In 


2  4 

particular,  for  j  =  2  and  j  =  4  we  must  have  y^  =  0  and  =  0, 


2  4 

respectively.  Consequently,  all  subsequent  subdivisions  of  ft  and  ft  must 


involve  an  index  different  from  I.  Otherwise,  an  optimal  solution  to  (5) 

J 


with  ft<=ftJ  (j=2,4)  would  have  both  x^.  and  y^  components  positive,  and  by 
Lemma  2  this  means  that  ft  is  fathomed.  Therefore,  any  path  in  the 


2  4 

branch-and-bound  tree  that  only  involves  subsets  of  the  form  0  and  ft 


can  have  at  most  n  branches. 


So  long  as  l  =  i  =  0,  the  index  i  can  conceivably  arise  repeat- 
Xi  yi 


edly  for  partitioning  a  sequence  of  nested  subsets.  We  argue  below  that 
the  sequence  is  finite  under  the  nondegeneracy  assumption.  Suppose  the 
index  I  is  used,  in  accordance  with  the  algorithm,  to  repeatedly  parti- 

i  ^  k  j 

tion  the  set  ftJ .  For  k=l,2,...,  let  ft  denote  the  subsets  of  ftJ  that 


4  4 

have  lower  bounds  of  zero  for  both  xT  and  yT;  that  is,  l  =  l  -  0. 

I  3 1’  XI  yi 


^  k  ^  k 

Similarly,  let  C  and  i|>  be  the  associated  sets  and  objective  functions 


in  problem  (5-j^).  Further  assume  that  a  noridegenerate  complementary 


*  *  ^  k 

solution,  say  (x  ,y  ),  is  in  C  for  all  k.  Thus,  we  assume  that  the 


*  * 


algorithm  is  converging  to  the  point  (x  , y  )  by  partitioning  ft  ,  for 


jk  Jk, 


each  k,  into  four  subsets  using  the  same  index  I.  With  (x  ,y  )  denot¬ 


ing  an  optimal  solution  to  problem  (5-j^)»  we  have 


1  ^1  ^2 
CJ  =  C  ^..., 


*  *  “ 

(x  ,y  )e  n  C  , 
k-1 


w-  v. 


IP 


* 

S*  *V  i 


kki 


m 


«'.  »-ST 
f  VNJ 


m 


m 


A 

V. 


V'V-I 

- 

'•V-W 


.  v, 


V  *  ‘ 
V  *0  4 

Wv,1 


m 


Jt  ^k.  *  * 

(x  ,y  )  +  (x  ,y  ), 


ik  jk 

where  we  assume  for  the  moment  that  the  sequence  {(x  ,y  )}  converges  to 

*  *  *  * 

(x  ,y  )  in  the  limit.  We  also  have,  by  construction,  that  (x  ,y  )  is  a 

Jk  *  * 

nondegenerate  extreme  point  of  SnC  for  every  k  whenever  (x  ,y  )  is  a 
nondegenerate  extreme  point  of  Snjj. 

*  * 

Under  the  nondegeneracy  assumption,  (x  ,y  )  is  a  locally  unique 
global  solution  to  the  problem 


minimize  x  y 


subject  to  (x,y)ESnC 


i  ^k  ^k  *  * 

for  every  k.  Since  {(x  ,y  )}  converges  to  (x  ,y  ),  then  for  suffi- 

*  * 

ciently  large  k  the  point  (x  ,y  )  uniquely  solves  the  linear  program 


.  .  ,  ,  JkNT  .  ,  ^k.T 

minimize  (y  ;  x  +  (x  )  y 


subject  to  (x,y)eSnC 


Thus,  convergence  is  achieved  finitely  when  a  linear  program  is  solved  at 

a  point  in  the  converging  sequence  which  is  sufficiently  close  to  the 

*  * 

nondegenerate  limit  point  (x  ,y  ).  This,  in  fact,  is  already  a  feature 
of  the  general  algorithm  in  [lj. 

The  following  scheme  is  implemented  in  the  algorithm  and  was  origi¬ 
nally  intended  to  accelerate  convergence  to  nonextreme  point  global  solu- 


r  v. 

/.y.y 


•W 

>  ■  •  * 


•  \ 

r/ 
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r.vwBST^y 


tlons  to  jointly  constrained  bilinear  programs.  Before  partitioning  ft 

jk  jk 

using  the  solution  (x  ,y  )  of  problem  (5-jk),  we  solve  the  linear  pro- 


Ji.  T  ^k  T 

minimize  (y  )  x  +  (x  )y 


subject  to  (x,y)eSnft 


k  T 

and  let  (x  ,y  )  denote  an  optimal  solution.  Finally,  we  minimize  x  y 

^  k  ^  k  k  k 

over  the  line  segment  connecting  (x  ,y  )  and  (x  ,y  ),  and  use  the 


minimizing  point  to  define  the  partiton  of  ft  .  Notice  that  the  only 

jk  jk 

difference  between  problems  (8)  and  (9)  is  in  the  sets  C  and  ft  which 


intersect  with  S.  In  general,  we  would  expect  the  optimal  objective 

jk  jk 

value  of  (9)  to  be  less  than  or  equal  to  that  of  (8)  because  C  cjj  . 

jk  jk 

However,  because  every  complementary  solution  in  ft  must  also  be  in  C  , 


for  k  sufficiently  large  and  assuming  nondegeneracy  of  the  limit  point 
*  *  ~^k  ^  k 

(x  ,y  ),  we  have  that  (x  ,y  )  also  solves  problem  (8).  Therefore,  our 
algorithm  converges  finitely  for  linear  complementarity  problems  that 


have  at  least  one  complementary  solution,  each  one  being  a  nondegenerate 
extreme  point  of  the  associated  polyhedral  set. 


3.  Concluding  Remarks 

A  related  enumeration  procedure  [3]  always  branches  into  two  nodes 


and  solves  the  subproblems  of  minimizing  x^  and  minimizing  y^  over  a  face 
(specified  by  the  parent  node)  of  the  set  S+  ■  {(x,y)cS:  x>0,  y>0),  where 


*  •  *  •  *  •  *  •  "  •  •  '  *  *  •  “  •  *  •  *  •  "  '  •  '  «  ’  •  *  «  •  »•  «  ■  » '  ■  •  .  •  •  «  *  .  •  .  *  ,  •  .  *  ,  •  •  •  »  '  "  •  '  .  *  .  '  _  -  ,  <  ,  •  _  *  . 
^  *^  ***•_?  *•*  l  '*•  »*.  ^ .  •**  •*.  f  •*  .  •*.  •*.  <*.  «* .  **.  *  •*,  •*.  ‘  V. 
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the  index  I  and  the  branching  node  are  selected  in  the  same  way  as  in  our 

algorithm  above.  If  the  minimum  value  of  x^  or  y^  is  zero,  then  the 

associated  node  represents  a  face  of  S+  having  a  lower  dimension  than  the 

face  of  the  parent  node;  otherwise,  the  corresponding  node  is  fathomed. 

A  comparison  between  the  bilinear  programming  approach  above  and  the 

enumeration  procedure  in  [3]  reveals  that  minimizing  x^  over  the  face 

4 

defined  by  the  branching  node  is  roughly  equivalent  to  minimizing  4/  (x,y) 
1  4 

over  S  n (C  uC  ).  Analogously,  minimizing  y  over  the  branching  node's 

2  12 

face  in  [3]  corresponds  to  minimizing  (x,y)  over  S  fi(C  u  C  )  in  our 
procedure.  Another  difference  between  the  two  strategies  is  that  the 
branching  variable  is  fixed  to  zero  in  descendant  nodes  (thus  defining 
the  faces  of  the  offspring)  in  [3],  whereas  the  branching  variable  in  our 
procedure  can  conceivably  take  on  nonzero  optimal  values  in  descendant 
subproblems,  although  this  is  discouraged  by  our  construction  which  main¬ 
tains  a  positive  objective  function  coefficient  for  the  branching  vari¬ 
ables  in  all  offspring  subproblems. 

An  implementation  of  the  bilinear  programming  procedure,  which  does 
not  exploit  the  refinements  described  herein,  was  tested  against  the 
enumeration  procedure  in  a  preliminary  experiment  reported  in  [2].  The 
enumeration  procedure  solved  all  seven  problems  faster,  but  only  two 
problems  ran  faster  by  one  order  of  magnitude  and  the  methods  had  run 
times  within  25%  of  each  other  on  another  two  problems.  This  suggests 
that  the  two  methods  are  comparable  on  some  problems. 

Despite  the  poor  implementation  of  the  bilinear  procedure,  it  never¬ 
theless  outperformed  the  cutting-plane  method  of  Ramarao  and  Shetty  [7] 
on  the  only  two  problems  in  common  among  the  test  problems  solved;  speci- 


fically,  problems  5  and  6  in  [2]  correspond  to  problems  PB-20-02  and 
PB-20-05  In  [7],  respectively.  These  problems  were  solved  on  the  same 
machine  (a  CDC  Cyber  170/730)  and  produced  run  times  favoring  our 
procedure  of  6.7  seconds  versus  27.57  seconds  for  problem  5  (PB-20-02) 
and  12.7  seconds  versus  25.85  seconds  for  problem  6  (PB-20-05).  Our  code 
was  unable  to  solve  problems  with  more  than  twenty  variables  because  of 
core  memory  limitations.  However,  the  memory  size  burden  can  be  con¬ 
siderably  reduced  by  the  refinements  to  the  procedure  presented  above. 

Another  novel  approach  to  solving  the  linear  complementarity  problem 
is  the  method  proposed  by  Solow  and  Sengupta  [8]  which  involves  finding 
roots  of  a  piecewise  linear  convex  function.  They  prove  their  algorithm 
to  be  finite  when  all  of  the  principal  minors  of  M  are  positive  (i.e., 
when  M  is  a  P-matrix).  While  their  approach  is  interesting,  the  computa¬ 
tional  results  they  reported  using  positive  definite  matrices  still  give 
Lemke's  method  [5]  an  edge.  By  contrast,  our  approach  is  finite  for 
arbitrary  matrices  and  can  solve  problems  on  which  Lemke's  method  fails. 
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